Astronomy & Astrophysics manuscript no. ggw03 


February 2, 2008 


(DOI: will be inserted by hand later) 





Monte Carlo radiative transfer 
in molecular cloud cores 

Jose Gongalves^'^, Daniele Galli^, and Malcolm Walmsley^ 

^ Centro de Astronomia e Astroflsica da Universidade de Lisboa, Tapada da Ajuda, 1349-018 Lisboa, Portugal 
■ ^ INAF-Osservatorio Astrofisico di Arcetri, Largo E. Fermi 5, 1-50125 Firenze, Italy 

o , 

, Received / Accepted 

(N . 

^ , Abstract. We present the results of a three-dimensional Monte Carlo radiative transfer code for starless molecular 

O ' cloud cores heated by an external isotropic or non-isotropic interstellar radiation field. The code computes the 

dust temperature distribution inside model clouds with specified but arbitrary density profiles. In particular 
we examine in detail spherical (Bonnor-Ebert) clouds, axisymmetric and non-axisymmetric toroids, and clouds 
heated by an external stellar source in addition to the general interstellar field. For these configurations, the code 
also computes maps of the emergent intensity at different wavelengths and arbitrary viewing angle, that can be 
^ . compared directly with continuum maps of prestellar cores. In the approximation where the dust temperature 

is independent of interactions with the gas and where the gas is heated both by collisions with dust grains 
\ and ionization by cosmic rays, the temperature distribution of the gas is also calculated. For cloud models with 

parameters typical of dense cores, the results show that the dust temperature decreases monotonically from a 
maximum value near the cloud's edge (14-15 K) to a minimum value at the cloud's center (6-7 K). Conversely, 
the gas temperature varies in a similar range, but, due to efficient dust-gas coupling in the inner regions and 
inefficient cosmic-ray heating in the outer regions, the gradient is non-monotonic and the gas temperature reaches 
a maximum value at intermediate radii. The emission computed for these models (at 350 /im and 1.3 mm) shows 
that deviations from spherical symmetry in the density and/or temperature distributions are generally reduced 
^ l' in the simulated intensity maps (even without beam convolution), especially at the longer wavelengths. 
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1. Introduction cal symmetry. More precisely, the maps of the millimeter 

dust emission show clear departures from circular sym- 

Understanding the structure of pre-protostellar cores is ^^^^^ ^^^^ ^ ^ ^^^^^ 2000; CaselH et al. 2002a,b; 

H ; an essential step towards an understandmg of protostellar ^^^^^^ 2002). In some cases there is also evidence 

, evolution. The density distribution immediately prior to polarization and hence for a magnetic field (Ward- 

the onset of gravitational collapse defines the initial con- Thompson et al. 2000). It seems likely to us that such 

ditions for collapse and hence one has a strong motivation behavior is caused by magnetic fields of energy densities 

to attempt to derive this density distribution from the sufficiently large to influence the core structure (see also 

observations. Since these objects are most commonly ob- ^gg^) ^^^^^ flattening along field lines (Basu 

served by means of their millimeter dust emission, one has 20OO, Jones & Basu 2002). We can attempt to infer the 

a strong interest m understanding the temperature distri- ^^^^-^^ distribution on the basis of mm-submm contin- 

bution of the dust m the pre-protostellar core and its mflu- ^^^^ ^^^-^^ ^^e temperature distribution 

ence on the emergent intensity distribution at millimeter- ^^^-^^^ ^^ere departures from spherical symmetry are 

submilhmeter wavelengths. In this article, we present a i^nportant. A first attempt in this direction was made by 

radiative transport code which has the interpretation of 2ucconi et al. (2001, hereafter ZWG, see also Evans et 

mm-submm maps of pre-protostellar cores as its mam 2OOI for the spherically symmetric case) but these au- 

S^^^- thors neglected heating of dust grains due to re-absorption 

The approach which we have adopted is influenced by photons emitted by the cloud itself. Here we present a 

the fact that, as a rule, the density structure of these Monte-Carlo code without this limitation but which is ca- 

objects shows evidence for large deviations from spheri- p^ble of handling a variety of geometries. We expect that 

this type of analysis will be especially useful as a tool for 
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the interpretation of the high class data we expect to come 
from future instruments like ALMA and HERSCHEL. 

Some high quality data arc already available and arc 
consistent with a picture in which the dust temperature of 
pre-protostellar cores decreases towards the center (Ward- 
Thompson ct al. 2002, Bianchi ct al. 2003). ISO obser- 
vations in the mid-IR (Bacmann et al. 2000) and in the 
far-IR (Ward- Thompson et al. 2002) suggest dust temper- 
ature gradients consistent with heating from the external 
interstellar radiation field (ISRF). The available models 
predict a factor of ~ 2 increase in temperature from cen- 
ter to edge (ZWG, Evans ct al. 2001, Andre ct al. 2003), 
with the gradient dependent on the cloud structure. The 
predicted emission of several cores is in reasonable accord 
with the data, and implies that the assumed equilibrium 
configurations are either unstable or maintained by a mag- 
netic field. 

Additional observational constraints on the thermal 
structure of prestellar cores are given by spatially resolved 
measurements of the gas temperature from NH3 observa- 
tions. Low-mass cores show fairly uniform gas temperature 
(e.g. Tafalla et al. 2002 for L1517B and L1498), whereas 
massive quiescent cores in Orion show significant temper- 
ature drops from edge to center (Li et al. 2003). Thus 
another check on our understanding of core structure is 
possible if radiative transfer models are also able to pre- 
dict the gas temperature distribution resulting from the 
balance of the relevant heating and cooling mechanisms. 
We therefore have also developed the capability to predict 
gas temperature distributions for our model cores. 

While techniques have been developed for the study 
of radiative transfer in cold cores in one dimension (e.g. 
Leung 1976, Rowan-Robinson 1980, Ivczic et al. 1997), 
their extension to three-dimensional geometries is not 
straightforward. A Monte Carlo technique has the advan- 
tage that it deals easily with any general geometry, inde- 
pendently of the density distribution and the anisotropy 
of the incident radiation field. The increasingly faster 
workstations available render this technique very useful, 
and fully 3-D Monte Carlo radiative transfer codes have 
been developed recently (e.g. Wolf et al. 1999, Niccolini 
et al. 2003). However, in these latter studies the radia- 
tive sources are internal, and this is a computationally 
distinct problem from a core heated from the outside for 
two main reasons. First, as we are interested in cold cores, 
it is reasonable to expect that they will be optically thin 
to their own emission, and therefore, for methods that 
use iteration, the temperature distribution quickly relaxes 
to the final value; alternatively, for methods that follow 
individual packets until they exit the domain, each indi- 
vidual packet will interact with the dust often only once. 
Second, as packets are launched from the outside, a wide 
range of spatial scales demands a prohibitively large pho- 
ton statistics, as the probability of a packet being launched 
in the direction of the innermost cells becomes increasingly 
small. 

In general, our paper confirms and supports the re- 
sults recently obtained for embedded prestellar cores by 



Stamatellos & Whitworth (2003) with an independent 
but similar Monte Carlo radiative transfer code. While 
Stamatellos & Whitworth (2003) focus on the important 
effect of the ambient medium in which (spherical) cores are 
embedded, our paper mostly emphasizes the consequences 
of deviations from spherical symmetry in the density dis- 
tribution and anisotropics of the ISRF. The two papers 
represent therefore complementary attempts to model ra- 
diative transfer in molecular cloud cores in less idealized 
situations than those considered so far. 

The structure of this paper is as follows. In Sect. 2, 
we set the general problem of radiative transfer in a dusty 
cloud. In Sect. 3, we describe the Monte Carlo method we 
have used. In Sect. 4, we present our choice of opacities and 
interstellar radiation field and briefly discuss their effect in 
the results. In esections 5 and 6 we present our results for 
two- and three-dimensional models, respectively. Finally, 
in Sect. 7 we summarize our conclusions. 

2. Radiative transfer in a dusty cloud 

In local thermodynamic equilibrium and ignoring the ef- 
fects of scattering of radiation by dust grains, the radiation 
field in a dusty cloud satisfies the equation of transfer 



n • V/.(r, n) = p(r)K,(r){B,[Td(r)] - J,(r, fi)}. 



(1) 



where I^, (r, n) is the specific intensity of radiation (in 
Jy ster~^) at position r in the direction fi at frequency 
u, p is the total (gas plus dust) density, Kj^ the total ab- 
sorption opacity, Td is the dust temperature, and Bi, is 
the Planck function. 



B.(Td) = 



2hu^ 



exp 
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The dust temperature distribution Td(r) is obtained by 
solving the equation of balance of emitted and absorbed 
radiation. 
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(4) 



is the mean intensity of radiation (in Jy). 

ZWG made the symplifying assumption that the dust 
grains are only heated by the incident interstellar ra- 
diation field Jl^^^ (assumed isotropic), neglecting re- 
emission of radiation absorbed by the grains themselves. 
That is to say, they assumed that the cloud is optically 
thin to its own radiation. Under this hypothesis, the equa- 
tion of radiative transfer has the solution 

jISRF r 

J^(r) = (p exp[-T^(r, n)] dn, (5) 



where 



pnR 



p(rr')/t^(fr') dr', 



(6) 
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is the optical depth from a point r inside the cloud to a 
point nR on the cloud boundary, in the direction fi. The 
calculation of the radiation field inside the cloud is thus 
reduced to a large number of integrations along different 
directions (rays) defined at each grid point. This simplifi- 
cation allowed ZWG to model axisymmetric non-spherical 
density distributions, but has the disadvantage of slightly 
underestimating the dust temperature, particularly at the 
cold center of the cloud. 

Our general approach is the following. First, for a given 
ISRF we compute the radiation field at any point inside 
the cloud core and solve the equation of balance of radi- 
ation emitted and absorbed by dust grains to obtain the 
dust temperature distribution rd(r). Second, we compute 
the flux emitted in a given view direction at different wave- 
lengths and we compare the results with observed spectral 
energy distributions and monochromatic maps for specific 
objects. 

In several cases we have also considered the effect of 
a spherical envelope surrounding the cloud core on the 
resulting temperature distribution. The rationale for this 
is that these objects are always embedded in a photon- 
dominated region (PDR) which transforms the incident 
optical-UV field into mid and far IR radiation. Often in 
practice, this is radiation from transiently heated parti- 
cles, either small grains or polycyclic aromatic hydrocar- 
bons, emitting in the 3-30 fim range. Since we wish to 
focus on the temperature structure of the high density 
core interior, we have decided to assume that the radiation 
field incident on embedded cores is cut off in the optical- 
UV wavelength range and that the radiation absorbed in 
the PDR is reradiated in the infrared. We thus incorpo- 
rate the photons reradiated by the PDR in the external 
field. The error involved here can be estimated by varying 
the form of the ISRF and appears to be small (though see 
the discussion of Andre et al. 2003). Our focus here is on 
the effects of the core geometry upon temperature struc- 
ture and we place minor emphasis on the spectral energy 
distribution of the incident field. 

We stress that our approach is general: the density 
distribution, the optical properties of the dust grains, the 
intensity and spectral energy of the radiation source(s) are 
input parameters that can be freely specified and can be 
adapted to model various astrophysical situations. 



3. The Monte Carlo method 

In the Monte Carlo technique, the computational domain 
is divided in a large number of cells of mass nii that absorb 
and emit radiation. The energy that enters the computa- 
tional domain is divided in N-y monochromatic packets of 
equal energy, that are launched stochastically from the 
boundary and followed until they exit the cloud. These 
packets may be absorbed by the dust particles, result- 
ing in a temperature increase of the cells where they are 
absorbed, and are immediately reemitted to enforce ra- 



diative equilibrium. If the cell i absorbs Ni packets, the 
energy absorbed by the cell per unit time is 



dE: 



abs 



dt 



\ rISRF 



(7) 



where L^^^^ is the luminosity of the ISRF at the cloud's 
surface, obtained integrating the mean intensity of the 
(anisotropic) ISRF over the cloud's surface and over fre- 
quency, 



(8) 



Assuming LTE conditions, the energy emitted by the i-th 
cell per unit time is 



dE? 



dt 



/>oo 

A-Kmi / K^B^{Td.i)dv 
Jo 



(9) 



where Td, i is the dust temperature of the i-th cell, obtained 
by equating the absorbed and emitted energies per unit 
time, 



/•OO 

/ Ki,By{TA 

JQ 



df 



-ISRF 



(10) 



The integral on the left-hand side of this equation can 
be tabulated for a grid of values of the dust temperature 
and opacity. The solution of Eq. (|10() can then be easily 
obtained by iteration and interpolation. In the limit to the 
continuum, — s- oo and — s- 0, Eq. (|10|l is equivalent 
to Eq. 0, but differs from the approximation adopted in 
ZWG, Eq. (jSJ, because in Ni are included the packets that 
may have been absorbed more than once in the cloud, and 
subsequently re-emitted. 

We handled re-radiation (and, therefore, energy con- 
servation) using the method devised by Bjorkman & 
Wood (2001): once a packet is absorbed, it is immedi- 
ately re-emitted (to conserve energy) at a new frequency 
determined by the local dust temperature. To achieve 
this, one starts by noting that if a cell has emissivity 
KiyB,y(Td,i — ATd,i) prior to absorbing a wave packet, then 
after packet absorption its temperature increases by an 
amount ATd.i, and the cell emissivity becomes k^B^ (Td,i). 
Thus, the increment in the cell emissivity is 



A> = K,[B,(Td,0-B.(rd,, - ATd,,)] 
'dB, 



dT, 



AT, 



d,i- 



(11) 



A photon packet is immediately reemitted with a fre- 
quency obtained from a probability distribution having 
the same spectral shape as Aj^, 



dPi _ Ky f dB^^ 
dv K \ dTd / rp^j' 



(12) 



where K — Ki,{dBy / dT^) dv is a normalization con- 
stant. The re-emitted photon packet is then followed until 
a new interaction occurs, and the procedure is repeated 
until the photon leaves the cloud. This method ensures 
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energy conservation in a statistical sense. Due to the low 
temperature of the cores, energy is generally re-radiated 
at long wavelengths, to which the core has a low optical 
depth, and the overall effect of re-emission on the dust 
temperature is more significant in the core's cold interior 
but is generally very small. 

3.1. Tests of the code 

We have tested our Monte Carlo code in different ways. 
Here we present the results of two complementary tests, 
one checking the accuracy of the packet propagation pro- 
cedure and one checking the equality of emitted and ab- 
sorbed energy. 

To test the packet propagation procedure, we have 
computed dust temperatures profiles for spherically sym- 
metric clouds, for which results can also be obtained 
with the ray integration method of ZWG. Specifically 
we have assumed the density profile of a Bonnor-Ebert 
(BE) sphere, an isothermal, pressure-bounded, equilib- 
rium model, that reproduces the centrally flattened den- 
sity distribution and the rapid density decrease at large 
radii observed in most cloud cores. Fig. ^ shows the dust 
temperature as function of radius for a BE with radius 
R = 0.1 pc, central density n(H2) = 4.4 x 10^ cm~^, and 
isothermal gas temperature Tg = 12 K (note that the gas 
temperature here is merely for the purpose of specifying 
the assumed density distribution). As shown in Fig.^ the 
result obtained using the ray integration method (solid 
line), is identical to the Monte Carlo result using the same 
approximation of optically thin re-emission (dotted line) . 
We also notice that when re-emission is taken into con- 
sideration (dashed line), there is only a small increase of 
the central temperature, even in this case, when the ex- 
tinction through the core center is ^2.2^m — 25, and the 
central dust temperature is little more than 6 K. 

To test energy conservation, we have computed the en- 
ergy emitted per unit time by an arbitrary area element 
on the cloud's surface by integrating I^, (the solution of 
eq. ^1 and /^^^^ over solid angles and frequencies, find- 
ing a very small discrepancy between the two values (less 
than 1 %). This test shows that for the given number and 
size of cells, the accuracy of Eq. (|1U|I and H12I) . that con- 
tain approximations introduced by the discretization of 
the computational volume, is very good. 

4. Dust opacities and the interstellar radiation 
field 

4.1. The interstellar radiation field 

Following ZWG and Evans et al. (2001), we adopt as a 
reference the ISRF given by Black (1994), which is an av- 
erage for the solar neighbourhood. Scaling the intensity of 
the ISRF by a factor Go leads to a change in the dust tem- 
perature at all radii by a factor of Gq^''^^^\ where /3 ~ 2 
is the power law index for the long wavelength opacity. 
Actually the radiation incident on molecular cloud cores 



may differ not only in intensity but also in spectral shape, 
depending on the degree of embedding of the cores in the 
ambient cloud and on the possible vicinity to hot stars. 
It appears for example to be the case that many cores 
associated with the p Oph star forming regions are sub- 
jected to an incident ISRF roughly an order of magnitude 
larger than the solar neighbourhood ISRF. An example of 
what can happen is given by a recent study by Andre et 
al. (2003) of a particular core, heated by a ISRF higher 
than the Black field by an order of magnitude in both the 
mid-IR and far-IR. This resulted in a higher central tem- 
perature and a sharper gradient at the edge than in our 
models for a similar density distribution. Thus observers 
attempting to infer density distributions from maps of the 
mm-submm dust emission must bear in mind the fact that 
both the radiation field and the geometry can play a role. 

In this paper however, our objective is to study the 
effects of geometry rather than those of the spectral 
characteristics or intensity of the incident radiation field. 
Geometry in this context can mean the geometry of the 
density distribution or that of the radiation field. We thus 
also consider a case where we use our Monte Carlo model 
to consider an anisotropic incident ISRF as exists in many 
galactic reflection nebulae. 

4.2. Dust opacities 

We adopt the dust opacities tabulated by Ossenkopf & 
Henning (1994, hereafter OH) in the wavelength range 
1 /im-1.3 mm. Recent work by Bianchi et al. (2003) and 
by Kramer et al. (2003) has shown that these theoretically 
derived opacities are consistent with the observed ratio of 
mm dust continuum intensity and near infrared extinc- 
tion. It is also clear from the study by OH that one can 
expect the dust opacity to vary between the lower den- 
sity outer layers of a core and the high density interior 
where most of the heavy element content of the core is in 
the form of ice of various sorts. We neglect such effects in 
this study mainly because, as shown below, the variation 
between different OH dust models is relatively mild. 

We label the opacity according to the column number 
in Table 1 of OH: OHl, standard MRN distribution; 0H4, 
MRN with thin ice mantles; 0H5, MRN with thin ice man- 
tles, after 10^ yr of coagulation at density 10^ cm""^; 0H8, 
MRN with thick ice mantles, after 10^ yr of coagulation 
at density 10^ cm~'^. 

The OH opacities are tabulated for wavelengths be- 
tween 1 pm and 1.3 mm. Within this range we have ob- 
tained values of the opacity at arbitrary wavelengths by 
four-point interpolation and above 1.3 mm, we have used 
a power law extrapolation based on the last two tabulated 
values. We neglect scattering completely as we are inter- 
ested in the transport of infrared photons. 

The effect of the various opacities on the temperature 
profile of a cloud is shown in Fig.|2 for the reference case of 
a BE sphere with central density n(H2) — 4.4 x 10^ cm^"^, 
radius 0.1 pc, gas temperature Tg = 12 K, embedded in a 
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spherical envelope with total extinction — 0.1. In 

this particular example, the gas-to-dust ratio was varied 
for each opacity prescription in order to keep fixed the 
total extinction at 2.2 fxia. 

As Fig. 121 shows, temperature variations due to differ- 
ent choices of the opacity are generally small. The opac- 
ities OHl, 0H4 and 0H5 have a similar dependence on 
wavelength, and the resulting temperature difference is 
equivalent to scaling the intensity of the ISRF by about 
10%. The 0H8 opacity, having a slightly different depen- 
dence on wavelength, results in a the temperature profile 
slightly different from the other three cases. In all remain- 
ing models presented in this study we have used 0H5. 

5. Two-dimensional models 

In this section we consider two particular situations that 
require a two-dimensional modeling of the radiative trans- 
fer. In the first example, the external radiation field is 
isotropic but the density distribution is not spherically 
symmetric. In the second example, the density profile is 
spherically symmetric but the incident radiation field is 
anisotropic, as the cloud is heated by an external stellar 
source. Both situations are easily handled by our Monte 
Carlo code. 

5.1. Singular isothermal toroids 

As a first application of our Monte Carlo code, we consider 
the radiative transfer in a singular isothermal toroid (Li & 
Shu 1996), a scale-free, axisymmetric equilibrium config- 
uration of an isothermal cloud under the influence of sef- 
gravity, gas pressure and magnetic forces. These toroids 
are characterized by a single parameter, Hq, representing 
the fractional amount of support provided by the magnetic 
field. We choose the particular value Hq = 0.5 to model a 
cloud with moderate axial ratio 2.5) intermediate be- 
tween a spherical unmagnetized cloud and a magnetically 
dominated disklike configuration. 

It should be noticed that, in all models presented, we 
are simply computing the dust temperature of pre-defined 
core density distributions, and therefore the resulting core 
is not the solution of the equilibrium equation if gas-dust 
coupling is assumed. Also, the practical definition of the 
boundary as an isobaric surface is only possible for spher- 
ical clouds, as in this case the isopycnic and isobaric sur- 
faces necessarily coincide. In all other geometries, we take 
the boundary as the isopycnic surface corresponding to a 
specified value of the density. 

In Fig. |31 are shown both the density distribution of 
the toroid and the resulting temperature distribution. The 
boundary of the configuration is defined as the isopycnic 
surface with n(H2) — 10* cm^'^. We choose this value 
of the density because it represents the threshold value 
that characterizes typical dense ("ammonia") cores. For 
a density profile characterized by a r"^ law, the visual 
extinction of the cloud material with density below this 
threshold value corresponds to about Ay — 2-3. Notice 



that the temperature varies over the boundary surface by 
a factor of about two from the origin to the maximum 
radius, and this difference is due to an increase of grain 
exposure to the ISRF as the radius increases, an effect 
known also for disks of low aspect ratio (e.g. Spagna et 
al. 1991), and is particularly enhanced in a toroid. It is 
also evident from Fig. |31 that the temperature decreases 
outward both in the horizontal and vertical directions. 

In Fig. 13 we also show synthetic (unconvolved with 
an observing beam) maps of the expected dust emission 
at 350 /im and 1.3 mm. In particular at 1.3 mm, the 
"observed contours" reflect fairly faithfully the projected 
column density and hence the density distribution. At 
350 /zm however, the intensity distribution is much more 
extended than at 1.3 mm. The half-power dimensions are 
~ 1.5" X 3.7" at 1.3 mm (in angular size at a distance 
of 140 parsec) but 18" x 25" at 350 fim. Thus the effect 
of heating from the exterior is mainly observable as an 
increase in size with frequency. 

5.2. Anisotropic Radiation Field 

While the above discussion concerns cores of non-spherical 
geometry subject to an isotropic incident radiation field, it 
should be noted that one expects often to find cases where 
the incident radiation field is anisotropic. Most obviously, 
this is the case in galactic reflection nebulae where the ra- 
diation field from single early type stars may be dominant. 
We therefore here (in a slight parenthesis to the discus- 
sion elsewhere in this paper) demonstrate the application 
of our code to such a case and compute the expected tem- 
perature distribution in a core subject to radiation from 
a nearby B star in addition to the isotropic ISRF (note 
earlier work in this field by Natta et al. 1981). 

The vicinity of a star to a prestellar core results in 
an anisotropic radiation field, leading to a 2-D radiation 
transfer problem in the case of a spherically symmetric 
core. It should be noted that this is a case where the ap- 
proximation of ZWG (Eq.EJ may not hold, as the expected 
increase in dust temperature may make the core no longer 
optically thin to its own radiation. 

To study the effect of such a field, we start by con- 
sidering a system composed of a prestellar core rep- 
resented as a BE sphere, with the physical parame- 
ters given in section 3.1, and a B3 star at a distance of 
0.15 pc from the edge of the core. The large range of spa- 
tial scales involved makes it difficult for a Monte Carlo 
program to solve this problem consistently, by emitting 
packets from the star, and following them through the 
ISM and the prestellar core. We avoid this problem by 
assuming that the stellar radiation is reprocessed by a 
tenuous but spatially extended photodissociation 
region (PDR) surrounding the core, and use a fit to 
the spectral emission of the PDR computed by Desert et 
al. (1990). In this way, the situation is reduced to the case 
treated in section 3.1, with the difference that the incident 
radiation is characterized by a PDR spectrum, instead of 
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the standard ISRF, and a variable intensity over the sur- 
face of the core. We further assume that the total lumi- 
nosity over any area element on the cloud's surface is the 
same as that obtained assuming no attenuation between 
the star and envelope. Thus, the luminosity of a surface 
element in spherical coordinates is 

dL(e) oc cos Oi sine d0, (13) 

where 9 is the angle between the point on the surface of 
the core and the star as seen from the core's center, r{9) is 
the distance of that point to the star, and 9i is the angle of 
incidence of the radiation. The proportionality constant is 
obtained by setting the total incident luminosity equal to 
(r2/47r)L*, where il is the solid angle of the core as seen 
from the star. 

The resulting temperature distribution is shown in the 
top panel of Fig. 0] As expected, the cloud is hotter on 
the side facing the star, and the temperature decreases 
both radially and azimuthally, with the lowest tempera- 
ture close to the centre of the core, which is almost twice 
higher than the minimum temperature of the Black ISRF 
heated core with the same density profile (Fig.^). The ra- 
tio L,/LiSRF _ 20, and the inverse square dependence 
of dL (see Eq. I13|l that makes this ratio over a surface ele- 
ment be much larger close to the star, has the consequence 
that the isotropic ISRF is negligible almost everywhere. 
The exception to this is on the core's side opposite to the 
star, where the shielding provided by the cloud's center 
suffices to make the isotropic ISRF the dominant heating 
field. 

In Fig. 21 we show the simulated emission maps at 
1.3 mm, 350 fim, 100 fim and 60 fim. The 1.3 mm map 
shows slight distortions due to the temperature structure 
but clearly allows a reasonable determination of the mass 
distribution in the BE sphere. As one goes to shorter wave- 
lengths however, the distortions become extreme and at 
60 microns, one sees essentially the heated surface of the 
cloud. Thus maps at different wavelengths may allow the 
identification of a heating source and show where the star 
is along the line of sight relative to the cloud. In such sit- 
uations, there will often be mid-IR data available which 
will additionally define the geometry of the "envelope" or 
PDR layer hypothesised in our analysis. 

6. Three-dimensional models 

Reverting to the case of an isotropic radiation field, we 
now consider the radiative transfer in a fully 3-D density 
distribution built on the results of Galli et al. (2001) for 
magnetized disks with uniform mass-to-flux ratio M (<&) /$ 
(isopedic) . A computationally convenient characteristic of 
these models is the existence of non-axisymmetric but an- 
alytical solutions of the equations of magnetostatics, cor- 
responding to surface density distributions S separable in 
polar coordinates {w^ip): 



where 



S(ro,(^) 



(15) 

1 + e cos if 

with < e < 1. Here, Q > 1 and e < 1 represent the 
increase of gas pressure and the reduction of the gravi- 
tational constant due to magnetic pressure and magnetic 
tension, respectively, related to the nondimensional mass- 
to-fiux ratio A = 2t:G'^/'^M{<^)/<^ by the relations (Shu & 
Li 1997) 



1 

A2' 



e 



A2 + 3 

A2 + r 



(16) 



We construct, in a rather ad hoc way, a 3-D model 
assuming a vertical stratification of gas in hydrostatic 
equilibrium and imposing the condition that the column 
density of the resulting density distribution reproduces 
Eq. ifT^ . obtaining 



S^icp) sech^ 



(17) 



where H{w, (p) — / iiGYilw , (p). 

As a specific example, we adopt the parameters de- 
rived by Galli et al. (2001) to match the thermal dust 
emission map obtained by Ward-Thompson et al. (1999) 
for the starless core L1544. We assume an effective sound 
speed a — 0.21 km s~^, a nondimensional mass-to-flux ra- 
tio A = 2, an eccentricity e = 0.54, and we define the core 
boundary by the isopycnic surface n(H2) = 10^ cm~'^. 

For this model we have also computed the gas temper- 
ature. To do this, we first notice that the energy deposited 
in the gas by cosmic ray ionization is negligible when com- 
pared to the energy absorbed by the dust, so that the en- 
ergy transfer between gas and dust will not significantly 
affect the grain temperature. Therefore, it is only neces- 
sary to solve the equation of thermal equilibrium of the 
gas. 



A,r - A„ 



0, 



(18) 



(14) 



with Fcr the cosmic-ray heating rate, Ag the gas cooling 
rate by molecular and atomic transitions, and Agd the gas- 
dust energy transfer rate. For these quantities, we assume 
the values given by Goldsmith (2001). We also assume 
a CO depletion factor of the form /dop = exp(n/ndop), 
where the critical density for CO depletion is taken to be 
fT-dop = 5.5 X 10*^ cm^'^. Fig. shows the density, dust 
temperature and gas temperature of our model, with the 
left and right panels showing the y = and the z = 
plane, respectively. 

This 3-D model and the 2-D toroidal model presented 
in the previous section show qualitatively similar temper- 
ature distributions, but the mismatch between isopycnic 
and dust isothermal surfaces is enhanced in the 3-D model, 
particularly in the y = plane where the geometry is more 
complex. 

The gas temperature, as it depends on both the dust 
temperature and on the gas density, shows a different dis- 
tribution from either (Fig. Ej). At the dense center, the 
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temperatures of the gas and the dust are well coupled, but 
they are decoupled at the lower density boundary. This 
results in a non-monotonic behaviour of the gas temper- 
ature, which reaches a maximum at intermediate density 
(an effect that also occurs in ID geometry, as shown in 
Galli et al. 2002). It is also evident, from Fig. (O that the 
gas temperature has a slightly higher value in the right 
side of the singularity than at the left, which is a con- 
sequence of the different dust temperature values at the 
density where the turnover in the gas temperature occurs 
on the two sides of the singularity. 

Fig. shows simulated dust emission maps at 1.3 mm 
and 350 /xm for the side view and top view. The flattening 
towards the center due to the dust temperature gradient is 
evident from the comparison between the emission at the 
two wavelengths. The vertical density stratification has a 
fairly strong effect on the shape of the map seen from the 
side, but is not noticable from the top. In general, one 
notes a slight "bottleneck" deviation from a purely ellip- 
tical shape of the singular isothermal disk model present 
in Galli et al. (2001). We note also that the "cometary 
shape" of the edge-on 1.3 mm map is rather similar to that 
observed towards L1544 and L63 by Ward-Thompson et 
al. (1999). 

7. Conclusions 

In this paper we have presented a Monte Carlo model for 
radiative transfer applicable to the study of the physical 
conditions in starless molecular cloud cores. The method 
described is able to compute accurately and efficiently the 
dust and gas temperature distribution in clouds of speci- 
fied but arbitrary density profiles heated by the ISRF, an 
external stellar source, or a combination thereof. In ad- 
dition, the method allows the computation of synthetic 
maps of the emission at near-infrared and submillimcter 
wavelengths for arbitrary viewing angles, and to predict 
the emergent spectral energy distribution. We have suc- 
cessfully tested our code by detailed comparison with one- 
dimensional calculations in spherical symmetry (ZWG, 
Galli et al. 2002) and have shown that the assumption of 
effective optical thinness is a good one for prestellar cores. 
As preliminary applications of our code, we have consid- 
ered the radiative transfer in model clouds of increasing 
degrees of asymmetry and in a spherical cloud core heated 
by a neighboring star. 

With our method, we confirm the results of ZWG, 
Evans et al. (2001), and Stamatellos & Whitworth (2003) 
that the dust temperature in cloud models with physical 
characteristics typical of dense cores varies monotonically 
between a minimum value at the center (6-7 K) and a 
maximum value near the cloud's edge (14-15 K) for the 
standard ISRF. 

In general, maps at submillimcter wavelengths closely 
reproduce the column density distribution, even for the 
most non-isothermal configurations. In the Ray leigh- Jeans 
regime, in fact, temperature differences of a factor ~ 2 
along the line-of-sight are washed out by the much larger 



variation of the density. However, the decrease of the dust 
temperature at the cloud center has the consequence that 
submillimcter maps sample preferentially the low-density 
envelope of prestellar cores, and are therefore unable to 
distinguish centrally peaked (pivotal) configurations from 
more flattened density profiles. We also note that our 
three-dimensional model based on the work of Galli et 
al. (2001) produces in natural fashion the cometary struc- 
ture observed towards many prestellar cores. 

For a spherical core heated by an external stellar source 
close to the cloud's surface, the temperature distribution 
is non-spherical, but the submillimcter emission generally 
reproduces the core's column density profile. Deviations 
from circular symmetry become extremely apparent on 
the other hand shortward of 100 /zm where one samples 
mainly the heated surface layer. 
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Fig. 1. Comparison of dust temperature profiles obtained 
with different methods for a BE sphere with central den- 
sity n(H2) = 4.4 X 10^ cm^'^, gas temperature Tg = 12 K, 
and radius R = 0.1 pc. The solid and dotted curves shows 
the result obtained with the ray integration method of 
ZWG and the Monte Carlo code, without including the 
effects of re-emission. The two results are in excellent 
agreement. The short- dashed curve shows the dust tem- 
perature computed with the Monte Carlo code including 
re-emission. In all three cases the cloud has no envelope. 
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Fig. 2. Dust temperature profiles (including re-emission) 
for the cloud model shown in Fig. 1 for different choices 

of the dust opacity: solid curve, OHl; dotted curve, 0H4: 
short-dashed curve, 0H5; long-dashed curve, 0H8 (see text 
for the definitions). The gas-to-dust ratio was varied so 
that the total extinction at 2.2 /im is the same in all cases. 
The cloud is embedded in a spherical outer envelope with 
total extinction A^J^jn =0.1. 
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Fig. 3. The top panel shows the dust temperature {top 
half) and density distribution {bottom half) of a singular 
isothermal toroid with Ha = 0.5, seen edge-on. The den- 
sity on the boundary is n{B.2) = 10^ cm~'^ and the isopy- 
cnic curves are logarithmically spaced by 0.5. The density 
distribution is azimuthally symmetric and also symmetric 
with respect to the midplane. The maximum outer tem- 
perature is 15 K and the isothermal curves are spaced by 
1 K. The bottom panel shows the emission at 1.3 mm {top 
half) and at 350 /um {bottom half). The isophotal curves 
are logarithmically spaced by 0.2 starting from the lower 
value 0.01 x I^'^^. With the adopted opacity law (0H5), 
the peak values are I^'^^ = 77 MJy sr~^ at 1.3 mm and 
382 MJy sr~^ at 350 /xm (in this and in the following fig- 
ures, the maps have not been convolved with an observing 
beam) . 
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Fig. 4. Model results for a spherical (Bonnor-Ebert) cloud 
core lieated by the standard ISRF and a B3 star lo- 
cated at 0.15 pc from the surface of the cloud (on the 
left side in this figure). The core has central density 
n(H2) = 4.4 X 10^ cm~^, gas temperature Tg = 10 K, and 
radius i? = 0.1 pc. The top panel shows the dust temper- 
ature {top half) and the density {bottom half). Isopycnic 
curves are logarithmically spaced by 0.5 and the bound- 
ary density is n(H2) = 2.3 x 10'^ cm~^. Isothermal curves 
are spaced by 2 K starting from = 14 K. The middle 
and bottom panels show the emission at 1.3 mm, 350 ^m, 
100 /im and 60 fxm. Isophotes are logarithmically spaced 
by 0.2 starting from 0.01 x J™^''. The peak intensity 7™'''' 
is 78 MJy sr"! at 1.3 mm, 3347 MJy sr'^ at 350 /xm, 
1404 MJy sr-i at 100 /xm and 200 MJy sr'^ at 60 /im. 
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Fig. 5. Density and temperature distribution in the magnetic 3-D model. Panels on the left (right) side show results 
in 2/ = (z = 0) plane, respectively. The top panels show the gas density, with isopycnic curves logarithmically spaced 
by 0.5, and density on the boundary n(H2) = 10* cm^'^. The middle panels show the dust temperature distribution, 
with isothermal curves spaced by 1 K between the minimum (7 K) and maximum (14 K) values. The bottom panels 
show the gas temperature, with isothermal levels ranging from 7 to 14 K. 
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Fig. 6. Normalized emission for the 3-D model, seen edge- 
on {top panel) and face-on {bottom panel). Each panel 
shows the emission at 1.3 mm {top half) and at 350 /xm 
{bottom half). Isophotal curves are logarithmically spaced 
by 0.2 starting from the lower value 0.01 x I™'^^. With the 
adopted opacity law (0H5), the peak values are: I^'^^ = 
89 MJy sr-i at 1.3 mm and 460 MJy sr'^ at 350 /im (edge- 
on case); I^^"" = 72 MJy st'^ at 1.3 mm and 307 MJy sr"! 
at 350 iJ.m (face-on case). 



